{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Planetary Ephemerides\n",
"\n",
"This tutorial demonstrates satkit's planetary ephemeris capabilities, covering both high-precision JPL ephemerides and faster low-precision analytical models.\n",
"\n",
"## Overview\n",
"\n",
"Satkit provides two approaches for computing solar system body positions:\n",
"\n",
"- **`sk.jplephem`** — High-precision positions from NASA/JPL Development Ephemerides (DE440/DE441). These are interpolated from pre-computed Chebyshev polynomial coefficients and deliver sub-kilometer accuracy across centuries.\n",
"- **`sk.sun` / `sk.moon`** — Low-precision analytical models from Vallado, useful when speed matters more than accuracy. The Sun model is accurate to ~0.01 degrees (1950--2050), and the Moon model to ~0.3 degrees in ecliptic longitude and ~1275 km in range.\n",
"\n",
"All positions are returned in meters in the GCRF (Geocentric Celestial Reference Frame) unless otherwise noted."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import satkit as sk\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"import scienceplots # noqa: F401\n",
"plt.style.use([\"science\", \"no-latex\", \"../satkit.mplstyle\"])\n",
"%config InlineBackend.figure_formats = ['svg']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Getting Planet Positions\n",
"\n",
"The `sk.jplephem` module provides geocentric and barycentric positions for all major solar system bodies. Bodies are specified using the `sk.solarsystem` enum.\n",
"\n",
"### Geocentric Positions\n",
"\n",
"`sk.jplephem.geocentric_pos` returns the position of a body relative to the Earth center, in the GCRF frame."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"t = sk.time(2026, 3, 20)\n",
"\n",
"# Sun position relative to Earth center\n",
"sun_pos = sk.jplephem.geocentric_pos(sk.solarsystem.Sun, t)\n",
"print(f\"Sun distance from Earth: {np.linalg.norm(sun_pos)/1e9:.3f} million km\")\n",
"print(f\"Sun position (GCRF): {sun_pos / 1e3} km\")\n",
"\n",
"# Moon position relative to Earth center\n",
"moon_pos = sk.jplephem.geocentric_pos(sk.solarsystem.Moon, t)\n",
"print(f\"\\nMoon distance from Earth: {np.linalg.norm(moon_pos)/1e3:.0f} km\")\n",
"print(f\"Moon position (GCRF): {moon_pos / 1e3} km\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Positions of other planets\n",
"planets = [\n",
" (\"Mercury\", sk.solarsystem.Mercury),\n",
" (\"Venus\", sk.solarsystem.Venus),\n",
" (\"Mars\", sk.solarsystem.Mars),\n",
" (\"Jupiter\", sk.solarsystem.Jupiter),\n",
" (\"Saturn\", sk.solarsystem.Saturn),\n",
"]\n",
"\n",
"for name, body in planets:\n",
" pos = sk.jplephem.geocentric_pos(body, t)\n",
" dist_au = np.linalg.norm(pos) / 1.496e11 # convert meters to AU\n",
" print(f\"{name:>8s}: {dist_au:.3f} AU from Earth\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Barycentric Positions\n",
"\n",
"`sk.jplephem.barycentric_pos` returns positions relative to the solar system barycenter. This is useful for interplanetary trajectory work or for understanding the geometry of the solar system.\n",
"\n",
"Note that the Sun's barycentric position is close to the origin but not exactly zero, since the barycenter shifts due to the gravitational influence of the giant planets."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Barycentric positions\n",
"sun_bary = sk.jplephem.barycentric_pos(sk.solarsystem.Sun, t)\n",
"earth_bary = sk.jplephem.barycentric_pos(sk.solarsystem.EMB, t)\n",
"jupiter_bary = sk.jplephem.barycentric_pos(sk.solarsystem.Jupiter, t)\n",
"\n",
"print(f\"Sun offset from barycenter: {np.linalg.norm(sun_bary)/1e3:.0f} km\")\n",
"print(f\"Earth-Moon barycenter: {np.linalg.norm(earth_bary)/1e11:.4f} x 10^8 km\")\n",
"print(f\"Jupiter: {np.linalg.norm(jupiter_bary)/1e11:.4f} x 10^8 km\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Vectorized Queries\n",
"\n",
"All `sk.jplephem` functions accept arrays of times, returning an Nx3 array of positions. This is much faster than calling the function in a Python loop."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create an array of times spanning one year\n",
"t0 = sk.time(2026, 1, 1)\n",
"times = np.array([t0 + sk.duration(days=d) for d in range(365)])\n",
"\n",
"# Get Mars geocentric position over the year (returns 365x3 array)\n",
"mars_pos = sk.jplephem.geocentric_pos(sk.solarsystem.Mars, times)\n",
"mars_dist = np.linalg.norm(mars_pos, axis=1) / 1.496e11 # AU\n",
"\n",
"print(f\"Mars position array shape: {mars_pos.shape}\")\n",
"print(f\"Mars distance range: {mars_dist.min():.2f} to {mars_dist.max():.2f} AU\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## State Vectors\n",
"\n",
"`sk.jplephem.geocentric_state` returns both position and velocity, which is essential for orbit determination and trajectory planning."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"t = sk.time(2026, 3, 20)\n",
"\n",
"# Moon state vector\n",
"pos, vel = sk.jplephem.geocentric_state(sk.solarsystem.Moon, t)\n",
"print(\"Moon state vector (GCRF):\")\n",
"print(f\" Position: [{pos[0]/1e3:.1f}, {pos[1]/1e3:.1f}, {pos[2]/1e3:.1f}] km\")\n",
"print(f\" Velocity: [{vel[0]:.3f}, {vel[1]:.3f}, {vel[2]:.3f}] m/s\")\n",
"print(f\" Speed: {np.linalg.norm(vel):.1f} m/s\")\n",
"\n",
"# Sun state vector\n",
"pos, vel = sk.jplephem.geocentric_state(sk.solarsystem.Sun, t)\n",
"print(\"\\nSun state vector (GCRF):\")\n",
"print(f\" Position: [{pos[0]/1e9:.4f}, {pos[1]/1e9:.4f}, {pos[2]/1e9:.4f}] x 10^6 km\")\n",
"print(f\" Velocity: [{vel[0]/1e3:.3f}, {vel[1]/1e3:.3f}, {vel[2]/1e3:.3f}] km/s\")\n",
"print(f\" Speed: {np.linalg.norm(vel)/1e3:.2f} km/s\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Low-Precision Models\n",
"\n",
"For applications where speed matters more than accuracy (e.g., quick visibility checks, rough force models), satkit provides analytical low-precision models:\n",
"\n",
"- **`sk.sun.pos_gcrf`** — Sun position accurate to ~0.01 degrees (Vallado Algorithm 29)\n",
"- **`sk.moon.pos_gcrf`** — Moon position accurate to ~0.3 degrees in ecliptic longitude and ~1275 km in range (Vallado Algorithm 31)\n",
"\n",
"These do not require JPL ephemeris data files and are faster to compute."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"t = sk.time(2026, 3, 20)\n",
"\n",
"# Low-precision Sun position\n",
"sun_lp = sk.sun.pos_gcrf(t)\n",
"print(f\"Sun distance (low-precision): {np.linalg.norm(sun_lp)/1e9:.3f} million km\")\n",
"\n",
"# Low-precision Moon position\n",
"moon_lp = sk.moon.pos_gcrf(t)\n",
"print(f\"Moon distance (low-precision): {np.linalg.norm(moon_lp)/1e3:.0f} km\")\n",
"\n",
"# Moon phase information\n",
"illum = sk.moon.illumination(t)\n",
"phase = sk.moon.phase_name(t)\n",
"print(f\"\\nMoon illumination: {illum*100:.1f}%\")\n",
"print(f\"Moon phase: {phase}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Comparing JPL vs Low-Precision Models\n",
"\n",
"How much accuracy do you sacrifice by using the analytical models? Let's compare Sun and Moon positions from both methods over an entire year and plot the position error."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Generate times over one year at daily intervals\n",
"t0 = sk.time(2026, 1, 1)\n",
"times = np.array([t0 + sk.duration(days=d) for d in range(365)])\n",
"days = np.arange(365)\n",
"\n",
"# JPL high-precision positions\n",
"sun_jpl = sk.jplephem.geocentric_pos(sk.solarsystem.Sun, times)\n",
"moon_jpl = sk.jplephem.geocentric_pos(sk.solarsystem.Moon, times)\n",
"\n",
"# Low-precision positions\n",
"sun_lp = sk.sun.pos_gcrf(times)\n",
"moon_lp = sk.moon.pos_gcrf(times)\n",
"\n",
"# Compute position errors\n",
"sun_err = np.linalg.norm(sun_jpl - sun_lp, axis=1) / 1e3 # km\n",
"moon_err = np.linalg.norm(moon_jpl - moon_lp, axis=1) / 1e3 # km\n",
"\n",
"print(f\"Sun position error — mean: {sun_err.mean():.0f} km, max: {sun_err.max():.0f} km\")\n",
"print(f\"Moon position error — mean: {moon_err.mean():.0f} km, max: {moon_err.max():.0f} km\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 7), sharex=True)\n",
"\n",
"ax1.plot(days, sun_err)\n",
"ax1.set_ylabel(\"Position Error (km)\")\n",
"ax1.set_title(\"Sun: Low-Precision vs JPL DE440\")\n",
"ax1.grid(True, alpha=0.3)\n",
"\n",
"ax2.plot(days, moon_err)\n",
"ax2.set_xlabel(\"Day of Year 2026\")\n",
"ax2.set_ylabel(\"Position Error (km)\")\n",
"ax2.set_title(\"Moon: Low-Precision vs JPL DE440\")\n",
"ax2.grid(True, alpha=0.3)\n",
"\n",
"fig.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Practical Example: Moon Distance Over a Month\n",
"\n",
"The Moon's orbit is noticeably elliptical, with perigee (closest approach) around 363,000 km and apogee (farthest point) around 405,000 km. Let's compute and visualize the Moon's distance over a month to identify perigee and apogee passages."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Compute Moon distance at hourly intervals over one month\n",
"t0 = sk.time(2026, 3, 1)\n",
"hours = np.arange(0, 30 * 24) # 30 days in hours\n",
"times = np.array([t0 + sk.duration(hours=float(h)) for h in hours])\n",
"\n",
"moon_pos = sk.jplephem.geocentric_pos(sk.solarsystem.Moon, times)\n",
"moon_dist = np.linalg.norm(moon_pos, axis=1) / 1e3 # km\n",
"\n",
"# Find perigee and apogee\n",
"perigee_idx = np.argmin(moon_dist)\n",
"apogee_idx = np.argmax(moon_dist)\n",
"\n",
"print(f\"Perigee: {moon_dist[perigee_idx]:.0f} km on {times[perigee_idx].datetime()}\")\n",
"print(f\"Apogee: {moon_dist[apogee_idx]:.0f} km on {times[apogee_idx].datetime()}\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"days_elapsed = hours / 24.0\n",
"\n",
"fig, ax = plt.subplots(figsize=(10, 5))\n",
"ax.plot(days_elapsed, moon_dist / 1e3, linewidth=1.5)\n",
"ax.plot(days_elapsed[perigee_idx], moon_dist[perigee_idx] / 1e3, 'v',\n",
" markersize=10, color='green', label=f'Perigee: {moon_dist[perigee_idx]:.0f} km')\n",
"ax.plot(days_elapsed[apogee_idx], moon_dist[apogee_idx] / 1e3, '^',\n",
" markersize=10, color='red', label=f'Apogee: {moon_dist[apogee_idx]:.0f} km')\n",
"\n",
"ax.set_xlabel(\"Days since March 1, 2026\")\n",
"ax.set_ylabel(\"Distance (x 1000 km)\")\n",
"ax.set_title(\"Earth-Moon Distance — March 2026\")\n",
"ax.legend()\n",
"ax.grid(True, alpha=0.3)\n",
"\n",
"fig.tight_layout()\n",
"plt.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"name": "python",
"version": "3.12.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}